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Abstract 

A new multiscale micromechanical approach is developed for the prediction of the behavior of fiber 
reinforced composites in presence of fiber clustering. The developed method is based on a coupled two- 
scale implementation of the High-Fidelity Generalized Method of Cells theory, wherein both the local and 
global scales are represented using this micromechanical method. Concentration tensors and effective 
constitutive equations are established on both scales and linked to establish the required coupling, thus 
providing the local fields throughout the composite as well as the global properties and effective 
nonlinear response. Two nondimensional parameters, in conjunction with actual composite micrographs, 
are used to characterize the clustering of fibers in the composite. Based on the predicted local fields, 
initial yield and damage envelopes are generated for various clustering parameters for a polymer matrix 
composite with both carbon and glass fibers. Nonlinear epoxy matrix behavior is also considered, with 
results in the form of effective nonlinear response curves, with varying fiber clustering and for two sets of 
nonlinear matrix parameters. 


Introduction 

Clustering in composite materials refers to the aggregation of the constituents such that the composite 
does not exhibit uniform microstructural distribution. Composite microstructures are sensitive to the 
specific manufacturing processes utilized, and because the processing cannot be perfectly controlled, 
some degree of clustering will always be present. Clustering is more prevalent in composites with low 
fiber or inclusion volume fractions as there is more volume available in which one constituent may be 
segregated. In fact, in nanocomposites, whose volume fractions are typically very low, clustering is a 
dominant phenomenon, which strongly influences the properties and performance of the composite 
material (c.f., Shaffer and Windle (1999), Vigolo et al. (2000), Shi et al. (2004)). 

Several investigators have presented methods to characterize the degree of clustering in composite 
materials (c.f., Guild and Summerscales (1993), Pyrz (1994), Scanlon et. al. (2003), Al-Ostaz et al. 

(2007), Vaughan and McCarthy (2010), Wilding and Fulwood (2011), Zangenberg et al. (2013)). In terms 
of modeling the effects of fiber clustering in composites, most investigations have been finite element 
based, as this approach can explicitly account for the composite microstructure. Examples include 
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Sorensen and Talreja (1993), Ghosh et al. (1997), Yang et al. (1997, 2000), Zeman and Sejnoha (2001), 
Wongsto and Li (2005), Melro et al. (2008), Abhilash et al. (2011), and Romanov et al. (2013). 
Additionally, the Generalized Method of Cells semi-analytical micromechanics theory has been employed 
to model random fiber distributions in continuous composites with a moving-window technique (Baxter 
and Graham (2000), Graham-Brady et al. (2003), Baxter et al. (2005), Acton and Graham-Brady (2010)). 

In the present investigation, the triply-periodic High-Fidelity Generalized Method of Cells (HFGMC) 
(Aboudi et al., 2013) is first enhanced to enable coupled multiscale analysis, wherein both the local and 
global scales are synergistically linked. The HFGMC method determines the strain (or stress) 
concentration tensors, which are used to establish the macroscopic constitutive equations of the 
composite, and also to provide the local stress and strain fields throughout the composite. This enables the 
prediction of not only effective composite properties, but also the effective nonlinear response (due to 
local inelasticity and damage) in response to full multiaxial loading, as well as generation of initial yield 
and damage envelopes. The present multiscale HFGMC implementation considers an arbitrary number of 
local repeating unit cells (RUCs), which consist of the composite constituent materials (e.g., fiber and 
matrix). These local RUCs are then arranged in a global RUC, whose behavior represents that of the 
multiscale composite material. The scales are explicitly coupled so that any local nonlinearity (which can 
be predicted using the local fields) is homogenized and passed to the global scale analysis. The examples 
presented here utilize a two scale analysis, but the methodology could be generalized to admit an arbitrary 
number of scales, as has been done in the case of the Multiscale Generalized Method of Cells (Liu et al., 
2011, Liu and Arnold, 2013; Bednarcyk et al., 2015). 

Herein, the multiscale HFGMC approach is employed to examine the effects of fiber clustering in 
continuous fiber reinforced composites. A key aspect of the approach taken involves segregating the 
composite microstructure into two zones, A and B, see Figure 1. Each zone consists of a composite 
material with its own fiber volume fraction, thus, by varying the fiber volume fractions and the 
arrangement/size of zones A and B, fiber clustering can be represented. For example, if the fiber volume 
fraction in zone A is higher than that in zone B, the fiber-rich zone A represents a clustered region. 
Conversely, if the fiber volume fraction in zone B is higher than that in zone A, zone A will be a matrix- 
rich region as compared to zone B. It should be noted that the present method admits an arbitrary number 
of zones. The alternative single scale HFGMC or finite element approach to modeling fiber clustering 
would involve full discretization of the actual composite microstructure. Obviously, this implies the use 
of a very dense mesh with great potential for meshing issues in the local fiber-matrix regions. In contrast, 
the proposed multiscale HFGMC methodology captures the effects of fiber clustering through local scale 
fiber volume fraction variations. Thus coarse meshes can be used on both scales, avoiding these meshing 
concerns while balancing efficiency and fidelity. Further, while mean field approaches (e.g., Mori-Tanaka 
method) can predict reliable effective properties, the predicted mean fields do not capture the variations of 
the stress fields in the composite constituents, even in the case of uniform fiber distribution. Utilization of 
mean fields to predict initial yield envelopes was shown to be problematic via comparison to finite 
element micromechanics analyses by Pindera and Aboudi (1988) in the case of uniform fiber distribution. 
In the presence of fiber clustering, the matrix stress field variations are clearly important, thus motivating 
the use of a higher fidelity micromechanical approach such as HFGMC. It should be noted, however, that 
lower fidelity methods with better computational efficiency, such as the Generalized Method of Cells or 
the Mori-Tanaka Method, could be used to generate the results like those presented herein. The impact on 
these results of the better local fields given by HFGMC versus other methods has not been investigated. 
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Figure 1. — Schematic drawing of a composite 
divided into two zones, A and B, which have 
different fiber volume fractions. 


In the present investigation, actual micrographs from polymer matrix composites are used to 
determine these two zones within the composite RUC. The fiber clustering is characterized by two 
nondimensional parameters, one representing the volume fraction of zone A within the composite, and the 
other representing the fraction of the fibers that are located in zone A. Studies are performed by varying 
these parameters and observing their impact on the predicted initial global-scale yield and damage 
envelopes, as well as the nonlinear global stress-strain response. It should be noted that the availability of 
the concentration tensors in HFGMC enables straightforward and efficient generation of these initial 
envelopes. In contrast, the use of the standard finite element method would require repeated application of 
many loading combinations to trace out these surfaces. 

The remainder of this paper is organized as follows. First, the multiscale HFGMC theory is presented 
for the generally triply-periodic case. Due to the possible nonlinearity of the constituent materials, an 
incremental (tangential) formulation is employed. The local and global concentration tensors and 
constitutive equations, as well as the coupling between them, are established. Next the methodology for 
determination of the initial yield and damage envelopes is described, followed by the method used to 
model fiber clustering. The Results and Discussion section details the processing of the composite 
micrographs to determine the analyzed composite RUCs. Results in the form of predicted effective 
properties, initial yield and damage envelopes, and macroscopic nonlinear stress-strain curves are 
presented for carbon/epoxy, glass/epoxy, and two different nonlinear representations of the epoxy matrix. 
The Conclusion section summarizes the paper and offers future research directions. 

Multiscale High-Fidelity Generalized Method of Cells 

In order to model the clustering of fibers in a composite material, the High-Fidelity Generalized 
Method of Cells (HFGMC) micromechanical model has been implemented in a two-scale framework. 
These two scales are referred to as global and local. The global scale represents a repeating unit cell of a 
periodic composite material, whose constituents are themselves periodic composite materials. Thus, the 
local scale represents the RUCs present within the global scale constituents, see Figure 2. 
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Figure 2. — (a) A multiphase composite with triply-periodic microstructure defined with respect to global coordinates 
(Xi, X 2 , X 3 ). (b) Its repeating unit cell (RUC) is represented with respect to the coordinates (Yi, Yz, Y 3 ) and is 
divided into A/a, A/b, and Nr global-scale subcells, in the Yi, Y 2 , and Y 3 directions, respectively, (c) The materials 
within the global-scale subcells (ABr) are represented by local-scale RUCs, defined with respect to local 
coordinates (yi, y 2 , / 3 ), and is divided into N a , A/p, and A/, local-scale subcells. 


Global Scale Analysis 

The HFGMC theory, which has been fully described by Aboudi et al. (2013), considers a composite 
material with triply-periodic microstructure, Figure 2(a), wherein periodicity conditions are enforced in 
all three Cartesian coordinate directions. The global repeating unit cell (RUC), Figure 2(b), defined with 
respect to local coordinates (Ti, Y 2 , TV), is divided into Na, Ms, and Nr global subcells in the Y\, Yi, and Y 3 
directions, respectively. Each global subcell is labeled by the indices (ABT) with A = 1 ,. . .,Na, 

B = 1,. . ,,Nb and T = 1,. . .,Nr, and may contain a distinct homogeneous material or a composite material. 
The dimensions of the RUC are D, FI, and L, whereas the dimensions of global subcell (ABF) in the Ti, 

Y 2 , and T 3 directions are denoted by Da, Ms, and Lr, respectively. A coordinate system ( T 2 ( A \ ) 

is introduced in each subcell whose origin is located at its center. The global subcell nonlinear elastic 
constitutive equation of the anisotropic material is given in an incremental form by, 
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where Acr| ABr) , Asj, ABr ^ , and C^ BI ) are the components of the stress increment, strain increment, and 
instantaneous (tangent) stiffness tensors of global subcell (ABT), respectively. 

The basic assumption in HFGMC is that the increment of the displacement vector AUj ABr) in each 

global subcell is represented as a second-order expansion in terms of its coordinates ( l) 1 A * , Y-j B 1 , l') 1 r * ) , as 
follows, 
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where Ae,y are the applied (external) average strain increments, and the unknown terms A must 

be determined from the fulfillment of the equilibrium conditions, the periodic boundary conditions, and 
the interfacial continuity conditions of displacements and tractions between global subcells. The periodic 
boundary conditions ensure that the increments of displacement and traction at opposite surfaces of the 
global RUC are identical. A principal ingredient in the HFGMC micromechanical analysis is that all these 
conditions are imposed in the average (integral) sense. 

As a result of the imposition of these conditions, a linear system of algebraic equations is obtained, 
which can be represented in the following form: 

K AV = AF (3) 

where the matrix K contains information on the global geometry and instantaneous properties of the 
materials (or composites) within the individual subcells (ABT), and the displacement vector increment, 

AV, contains the unknown displacement coefficients AJF^^p , which appear on the right-hand side of 
Equation (2). The vector AF contains information on the applied average strain increments As )j . The 
solution of Equation (3) enables the establishment of the following localization relation which expresses 
the average strain increments As^ 1 * 11 in the global subcell (ABT) to the externally applied average 
strain increments As ;/ in the form, 


< ,n =4“ n AS* (4) 

where ) arc the instantaneous global strain concentration tensor components, of the subcell (ABT). 

The final form of the effective incremental constitutive law of the multiphase composite, which 
relates the average stress increments Act ;/ - and strain increments As , is established as follows: 
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In this equation are components of the instantaneous effective global stiffness tensor, which are 
given by, 
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Next, the components of the global instantaneous stress concentration tensor, , which relate the 

average stress increments in the global subcell, Aa-y ABr - ) , to the average (global) stress increments, Ao)y , 
are determined. By combining Equations (1) and (4), the subcell stresses are given by, 
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and S* are components of the effective instantaneous global compliance tensor, which are the 
components of the inverse of the instantaneous global effective stiffness tensor, C*, see Equation (6). 


Local Scale Analysis 

As stated above and shown in Figure 2, the global subcells can be occupied by periodic composite 
materials, whose responses are determined by the micromechanical analysis of local RUCs. Therefore the 
global scale analysis subcell quantities are the RUC quantities in the local scale analysis. As such, the 
global subcell incremental constitutive equation (1) serves as the RUC constitutive equation for the local 
scale. The local subcell incremental nonlinear elastic constitutive equation is given by, 


A(j (ABr,ap Y ) = c (ABr,afr) ^(ABEaPy) 
y ijkl kl 


( 10 ) 


where the indices (ABT) indicate the global subcell and the indices (aPy) indicate the local subcell. As in 
the global scale analysis, the increment of the displacement vector A« ; - ABr,a ^ y) in each local subcell is 

expanded into quadratic forms in terms the coordinates ( A B r ’ u ’ , y B r 19 ^ jt* A B r A ) j centered in the 
local subcell, as follows, 
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where p = A,B,r for j = 1,2,3, respectively, and d( ABr,a), /a abi ,p>, and /(ABr,y>are the dimensions of the local 
subcell (aPy) within the global subcell (ABT). The unknown local terms are determined 

from the fulfillment of the equilibrium conditions, the periodic boundary conditions, and the interfacial 
continuity conditions of displacements and tractions between local subcells. Again, all of these conditions 
are imposed in the average sense. 

The application of these conditions lead to the local system of linear algebraic equations, which can 
be written in the following form: 


k Av = Af (12) 

where the matrix k contains information on the local geometry and instantaneous properties of the 
materials within the individual subcells (aPy), and the local displacement vector increment, Av, contains 
the unknown local displacement coefficients Aw^^y 01 ^ ^ which appear on the right-hand side of 
Equation (11). The vector Af contains information on the global subcell average strain increments 
Ae^ ABr * . The solution of Equation (12) enables the establishment of the following localization relation 

which expresses the average local strain increments As( ABr,a ^ / * in the local subcell (aPy) to the global 
subcell average strain increments As^ ABr> in the global subcell (ABT), 

A -(ABr,aPy) =J (ABr,aPy) A -(ABr) (13) 


where are the instantaneous strain concentration tensor components of the local subcell (aPy) 

within the global subcell (ABT). With the established ^4(^ Br,(X ^ , the local RUC effective instantaneous 

stiffness tensor, which corresponds to the global subcell effective instantaneous stiffness tensor, c( ABr - ) , 
in Equation (1), can be written as, 
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As before, the components of the local instantaneous stress concentration tensor, , which relate 

the average stress increments in the local subcell, Aaj / ABr,a ^, to the global subcell average stress 
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increments, Aa ? y ABr ^ , are determined. By combining Equations (10) and (13), the local subcell stresses 
are given by, 
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and SIT' are components of the effective instantaneous global subcell compliance tensor, which 

are the components of the inverse of the instantaneous global subcell stiffness tensor, C (ABr) , see 
Equation (14). 

Therefore, starting at the local scale, given the properties of the local-scale constituents (e.g., fiber 
and matrix) and their arrangement in the RUC, the local HFGMC provides the incremental system of 
algebraic equations (12), from which the instantaneous strain concentration tensor, A (Alil - ttf:,7) , is 
determined. Equation (14) then provides the instantaneous effective stiffness tensor, C (ABr) , of each 
composite material on the local scale . At the global scale, these materials are located within the global 
subcells (ABT). Given the arrangement of these composite materials within the global RUC, the global 
HFGMC analysis provides the incremental system of algebraic equations (3), from which the global 
instantaneous strain concentration tensor, A (ABr) , is determined. Then, the effective instantaneous stiffness 
tensor of the entire multiscale material, C*, can be determined from Equation (6). Given external applied 
loading at the global scale, the stress and strain fields throughout the multiscale material (at both scales) 
can be obtained through the localization equations, Equations (4) and (8) for the global scale and 
Equations (13) and (16) for the local scale. 


Initial Yield and Damage Envelopes — Methodology 

In the present section, the methodology used to generate initial yield and damage envelopes within 
the multiscale HFGMC analysis is described. The yield surfaces are based on the traditional von Mises 
criteria, as defined by the equivalent stress in the local subcell, given by, 


-(ABr,«Py) _ 

where a jy ABr ' u ^ ; ' ) = q! AB r > a Pt > - a| ABr ’ tt l ! ' ; ' ) ( 5 // /?> are the subcell deviatoric stresses at the local scale, 5,, 

is the Kronecker delta, and Y is the yield stress of the material in the local subcell in simple tension. 

The damage initiation criterion, suggested by Lemaitre and Chaboche (1990), is used to express the 
critical strain energy release rate associated with loss of stiffness in brittle materials, which are typically 
highly dependent on the state of triaxial stress, a/, = Gm 73 . The average triaxiality function (Lemaitre, 
2001) in the local subcell is given by, 


3_£(ABl\apy)£(ABr,aPy) _ 

2 *7 7 


= 7 


(18) 


NASA/TM— 2015-218838 



( 19 ) 


n (ABr.aPy) 

/v v 


2 

I 


f 

(1 + v ( ABr ’«Py) )+ 3(1 - 2vi ABr >«Py) ) ■ 

v 


-(ABr.aPy) ' 

— (ABr,aPy) 
°eq 


2 


where v (Alil uf:,7 > is the Poisson’s ratio of the isotropic local subcell and CT AABr ’ a P'/) j s the local subcell 
average hydrostatic pressure. The damage initiation criterion is then given by, 


-(ABr,apy) /w (ABr.apy) 

®eq v ® cr 


( 20 ) 


Where a c , is a material parameter that specifies the stress at which damage initiates for a uniaxial tensile 
test. 

To generate a yield or damage initiation envelope in a given global scale stress plane, for example, 
cfi | - a 22 , the global stress component increments are expressed as, 


Aajj^sina AG ?2 =^'cosa (21) 

where R is the radial distance from the origin to a point located on the initial yield or damage surface in the 
cti 1 -a 22 stress plane and a is the corresponding polar angle. For a given a, the global subcell stress 

increments, Acrj ABr ^ , are readily determined (in terms of R) by substituting Equation (21) into Equation (8), 


Aaf Br > = 4 ABr) R sina + B ^ R cosa 


( 22 ) 


Substituting this equation into the local stress concentration equation ( 1 6), 

^ABr*, =M ( r ,„ Pr) ( 5 (ABn sma +5 <ABr, coso ) (23) 

Due to initial linearity, the yield or damage initiation envelope can be reached in just one increment, 
therefore, the local subcell stresses can be considered as equivalent to their increments. Then the local 
subcell stresses given by Equation (23), which are proportional to R, can be utilized in either initiation 
criterion. Equation (18) or (20), thus providing the value of RJY or R/<5„- corresponding to yield or damage 
initiation at the specified polar angle a. 


Modeling of Fiber Clustering 

Herein, the fiber clustering is modeled using two zones with two different volume fractions, which 
may be arranged arbitrarily. To this end, let the fiber reinforced composite material be composed of two 
distinct regions, zone A and zone B, which are distributed within the composite in an arbitrary manner, 
see Figure 1. Note that the number of zones is arbitrary, but herein we have chosen to use two to 
maximize the effect of clustering. Following Shi et al. (2004), the volume fractions of the fibers within 

each zone are denoted as cj and c'j . Let V be the total volume of the composite and V A the volume of 

zone A. In addition, let Vf and V 4 be the total volume of fibers in the composite and within zone A, 

respectively. Consequently, the overall (average) volume fraction of the fibers in the composite is, 
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Further, the volume fraction of zone A within the composite is given by, 


V A 



and the fraction of the fibers that are located in zone A is given by, 
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These equations imply that within zone A, the fiber volume fraction is, 


( 25 ) 
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and within zone B, the fiber volume fraction is, 



1 ~c A 


(27) 


(28) 


Note that, when ca = C, the fibers are uniformly distributed in the composite, and thus no clustering is 
present. When ca = 1 , the entire composite is occupied by zone A, whereas, when ca = 0, the entire 
composite is occupied by zone B, with both of these conditions also representing no clustering. When 
C, = 1 , all fibers are located in zone A, while zone B is pure matrix, which represents complete clustering. 
Likewise, when C, = 0, all fibers are located in zone B, while zone A is pure matrix, which also represents 
complete clustering. Therefore, for a given segregation of the composite into zones A and B, the degree of 
fiber clustering can be fully described by the values of ca and C, both of which vary between 0 and 1 . 

For a given composite overall volume fraction, c/, and a given volume fraction of zone A within the 
composite, c A , the above equations impose limits on the variation C. In order for the fiber volume fractions 
of zones A and B, c4 and , to remain less than 1 , the maximum and minimum values of C are given by, 


Cmin — 


Cf+C A - 1 

C f 


_CA_ 
max — 

c r 


(29) 


Results and Discussion 

The effect of fiber clustering could be investigated by arbitrarily selecting the global scale RUC and 
choosing various values of ca and C. Flere, on the other hand, a micrograph of the 0° ply within a continuous 
fiber carbon/epoxy composite, see Figure 3, has been used to provide a more realistic global scale RUC. 

The overall fiber volume fraction in this micrograph, determined via image pixel analysis is c/= 0.457. To 
assist in dividing the composite into zones A and B, a 10 by 16 square grid was placed over the micrograph, 
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and the fiber volume fraction of each square was determined via image pixel analysis, see Table 1. Then, the 
grid squares could be distributed to zones A and B based on a cutoff value of the fiber volume fraction. That 
is, given a cutoff value for the fiber volume fraction (e.g., 0.55), all grid squares above this value were 
assigned to zone A, while the remaining squares were assigned to zone B. While the choice of cutoff value 
is arbitrary, herein three such values (0.55, 0.45, 0.38) were chosen to investigate the effect of fiber 
clustering, see Figure 4. Note that, once a cutoff value is chosen, the value of ca is determined as the fraction 
of shaded (zone A) squares in Figure 4. The corresponding values of ca are 0.26, 0.5, and 0.74, respectively. 
The value of C, on the other hand, must be specified as constrained by Equation (29). Thus, for ca = 0.26, 

0 < C, < 0.56, for ca = 0.5, 0 < C, < 1, and for ca = 0.74, 0.44 < C, < 1. Once C, was chosen, the fiber volume 
fractions within zones A and B were calculated from Equations (27) and (28). 



20 pm 


Figure 3. — Micrograph from the 0° ply in an IM7/MTM-45 
carbon/epoxy composite. The overall fiber volume fraction, 
based on image pixel analysis, is 0.457. 


TABLE 1.— FIBER VOLUME FRACTIONS OF THE 10 BY 16 SQUARE GRID THAT DIVIDES THE MICROGRAPH OF 
FIGURE 3. THE SHADED ENTRIES ARE THOSE WITH FIBER VOLUME FRACTIONS ABOVE 0.55, WHICH 
CORRESPOND TO THOSE SHADED IN FIGURE 4(A) 



i 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

1 

0.648 

0.507 

0.715 

0.701 

0.609 

0.466 

0.562 

0.485 

0.573 

0.449 

0.411 

0.431 

0.328 

0.577 

0.271 

0.267 

2 

0.509 

0.512 

0.435 

0.506 

0.577 

0.575 

0.706 

0.604 

0.662 

0.714 

0.592 

0.645 

0.638 

0.650 

0.382 

0.417 

3 

0.447 

0.473 

0.357 

0.465 

0.518 

0.483 

0.559 

0.469 

0.548 

0.537 

0.447 

0.387 

0.428 

0.482 

0.450 

0.446 

4 

0.380 

0.490 

0.554 

0.252 

0.487 

0.381 

0.537 

0.611 

0.586 

0.532 

0.404 

0.362 

0.516 

0.314 

0.373 

0.454 

5 

0.375 

0.515 

0.612 

0.627 

0.674 

0.602 

0.559 

0.573 

0.573 

0.573 

0.652 

0.468 

0.486 

0.305 

0.362 

0.271 

6 

0.278 

0.384 

0.366 

0.404 

0.515 

0.692 

0.699 

0.652 

0.654 

0.566 

0.570 

0.642 

0.424 

0.310 

0.234 

0.299 

7 

0.188 

0.182 

0.253 

0.110 

0.419 

0.438 

0.347 

0.359 

0.421 

0.598 

0.489 

0.409 

0.306 

0.342 

0.200 

0.317 

8 

0.068 

0.341 

0.507 

0.373 

0.294 

0.507 

0.395 

0.242 

0.359 

0.341 

0.385 

0.386 

0.436 

0.369 

0.464 

0.506 

9 

0.298 

0.414 

0.382 

0.386 

0.385 

0.465 

0.596 

0.520 

0.455 

0.565 

0.462 

0.374 

0.361 

0.440 

0.406 

0.385 

m 

0.446 

0.482 

0.389 

0.410 

0.418 

0.442 

0.416 

0.472 

0.451 

0.439 

0.330 

0.500 

0.330 

0.534 

0.597 

0.462 
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□ zone A 

□ zone B 


Figure 4. — Division of the micrograph (Figure 3) into zones A and B based on a cutoff value of the fiber volume 
fraction, (a) Cutoff fiber volume fraction = 0.55 and ca = 0.26, (b) Cutoff fiber volume fraction = 0.45 and ca = 0.5 
(c) Cutoff fiber volume fraction = 0.38 and ca = 0.74. 


The microstructures shown in Figure 4 were employed as the global RUCs in the multiscale HFGMC 
analyses. At the local scale, 3 by 3 RUCs were employed to maximize efficiency, where the fiber volume 
fraction in zones A and B were adjusted to obtain the desired local fiber volume fractions, see Figure 5, 
Figure 6, and Figure 7. Note that square fiber packing was assumed at the local scale throughout this 
investigation. Table 2 provides a summary of the analysis cases considered associated with the 
micrograph in Figure 3. 

The transversely isotropic carbon fiber properties used within the local scale RUCs are as 
follows: axial Young’s modulus = 286 GPa, transverse Young's modulus = 12.4 GPa, axial shear 
modulus = 20 GPa, axial Poisson’s ratio = 0.29, and transverse Poisson's ratio = 0.29. The isotropic epoxy 
matrix properties considered are: Young’s modulus = 5.2 GPa and Poisson's ratio = 0.35. Throughout this 
study, the continuous fibers are oriented in the one-direction. 
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□ Fiber 
] Matrix 


c A f = 0.89 


c B f — 0.31 


-...x c?=0.18 c/= 0.55 

Figure 5. — Multiscale material distribution for ca = 0.26 with (a) C, = 0.5 and (b) C, - 0.1. 




Fiber 

Matrix 


0.68 c B = 0.24 cj= 0.22 c}= 0.70 

Figure 6. — Multiscale material distribution for ca = 0.5 with (a) C, = 0.74 and (b) C, = 0.24. 
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□ Fiber 
] Matrix 


cj = 0.55 7=0.20 


cj = 0.30 0.93 


Figure 7. — Multiscale material distribution for ca = 0.74 with (a) C, = 0.89 and (b) C, = 0.48. 


TABLE 2.— SUMMARY OF THE FIBER CLUSTERING PARAMETER CASES CONSIDERED 
FOR THE COMPOSITE MICROGRAPH SHOWN IN FIGURE 3 


Cut off fiber volume 
fraction 

CA 


C A 

c f 

C B 

c f 

Used to divide squares in 
micrograph into 
zones A and B 

Volume fraction of 
zone A within the 
composite 

Fraction of the fibers 
that are located in 
zone A 

Fiber volume fraction 
within zone A 

Fiber volume fraction 
within zone B 

0.55 

0.26 

0.5 

0.89 

0.31 



0.1 

0.18 

0.55 

0.45 

0.5 

0.74 

0.68 

0.24 



0.24 

0.22 

0.7 

0.38 

0.74 

0.89 

0.55 

0.2 



0.48 

0.3 

0.93 


Effective Elastic Properties 

The effect of fiber clustering on the effective properties can be extracted from the effective stiffness 
tensor, C*. In fact, the determination of C* does not require the coupled multiscale approach employed 
herein. Rather, by homogenizing zones A and B independently, and then employing the resulting stiffness 
tensors in a separate, uncoupled, homogenization procedure, C could also be determined. Here, however, 
the fully coupled approach was used. It is emphasized that the initial yield and damage envelopes 
presented in the next section, based on constituent properties, necessitate the use of a coupled multiscale 
approach as the fields at the local scale are required. 

Figure 5 shows the microstructural details of the case where ca = 0.26. This figure illustrates the 
global and local RUC geometries used to consider the cases of C, = 0.5 and C, = 0.1. For C, = 0.5, one half 
of the fibers are located in zone A, which is considerably smaller than zone B. Therefore, the volume 
fraction in zone A is quite high ( cj = 0.89), whereas cj = 0.31 (maintaining the overall c/= 0.457). In 
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contrast, for C = 0.1, only 10 percent of the fibers are located in zone A. Thus the fiber volume fractions 
of zone A and zone B are reversed, with cj =0.18, whereas = 0.55. Both cases represent deviations 

from uniform distribution of fibers in the composite; i.e., clustering. The additional cases where ca = 0.5 
and 0.74 are depicted in Figure 6 and Figure 7. 

Figure 8 shows the normalized effective axial and transverse (in both directions) Young's moduli (E n, 
E 22 , and £ 33 ) and the normalized effective axial shear modulus ( G 1 2 ) of the composite for different values 
of ca and £. For each of the three cases in Figure 4, two values of the parameter C, which indicates the 
fraction of fibers located in zone A, are considered: one which provides higher fiber volume fraction in 
zone A and the other that provides higher fiber volume fraction in zone B. The results are normalized by 
the uniform fiber distribution case, which can be obtained by setting ca = C- Clearly, the impact of fiber 
clustering on the effective properties of the composite is minor. The maximum variation, which occurred 
for Gn, is approximately 6 percent versus the uniform fiber distribution case. This validates the traditional 
approach of predicting composite effective elastic properties by considering a uniform fiber distribution. 

It is interesting to note that, while clustering resulted in slight reductions in the effective transverse 
Y oung's moduli, E 22 and £33, it resulted in a slight increase in the effective axial shear modulus, Gn. Of 
course, the comparison is made to the case of uniform fiber distribution, which herein assumes square 
fiber packing. 


Initial Yield and Damage Envelopes 

In this section, the initial yield and damage envelopes are presented for several global scale stress 
planes while varying the fiber clustering. A parametric study is performed by varying the value of ca as 
shown in Figure 4. As in the case of the effective elastic properties above, for each of the three cases in 
Figure 4, two values of the parameter C, are considered: one which provides higher fiber volume fraction 
in zone A and the other that provides higher fiber volume fraction in zone B. Details are shown in 
Figure 5 to Figure 7. Also shown in the initial envelopes for comparison is corresponding case with 
uniform fiber distribution. 

1.08 -| 

C A = 0.26 

• C = 0.1 

□ <^=0.74 

C A = 0.5 
a ^=0.24 

□ ^ =0.74 

C A = 0.74 

□ 0.48 

□ <^= 0.89 


Ell E22 E33 G12 

Figure 8. — The effect of fiber clustering on the normalized effective elastic properties of the composite. 
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Figure 9 shows the initial yield and damage envelopes for ca = 0.26 for uniform fiber distribution 
and clustering with C, = 0.5 and C = 0. 1 . Note that the stress magnitudes are normalized with respect to 
the yield stress ( Y m ) or critical damage stress ( a cr -m ) of the isotropic matrix in simple tension, see 
Equations (18) and (20). Further, the predictions of the initial envelopes are based on the yielding and 
damage of the first matrix subcell within the local scale RUCs. That is, the fiber is not permitted to yield 
or damage. 

Comparing the initial yield envelopes with their corresponding initial damage envelopes, the effect of 
the inclusion of hydrostatic stress in the damage criterion is seen as the damage envelopes are in all cases 
smaller than the corresponding yield envelopes. This is particularly true in the 022-033 envelopes, wherein 
the global hydrostatic stress is maximized in quadrants 1 and 111. Clearly, the effect of clustering is 
significant since the envelopes that include this effect are always inside the envelopes associated with 
uniform fiber distribution. It should be emphasized that these initiation envelopes may or may not be 
representative of composite ultimate strength. In brittle systems (e.g., ceramic matrix composites), 
however, they would be expected to be more representative of final failure. In addition, the C, = 0.5 
envelopes are always inside the C = 0. 1 envelopes. This is likely caused by the very high fiber volume 
fraction of zone A {c4 = 0.89) in the case of C, = 0.5, which results in high local scale stress 

concentrations in the matrix. 

It is interesting that initiation of yielding and damage may occur in either zone A or zone B at each 
point on each envelope, thus diamond symbols are included in Figure 9 to indicate the portions of each 
envelope where initiation occurred in zone A. Transitions between zone A initiation and zone B initiation 
occur in many of the presented envelopes, while some exhibit initiation only in zone B (solid lines with 
no symbols). Obviously, for the case of uniform fiber distribution, there is no distinction between zones A 
and B, thus such transitions are not applicable. 

Examining the initial 01 1-022 yield and damage envelopes reveals that, as expected, the effect of the 
continuous fibers is dominant as there is approximately a 25:1 difference in initiation stress for pure an 
versus pure G 22 . Also, the effect of clustering is negligible under pure axial loading. Also of note is the 
fact that the damage initiation always occurs in zone B irrespective of the value of C,. In contrast, this 
parameter has a pronounced effect on the initiation location in the case of yielding. 

Considering biaxial loading in the transverse directions ( 033 - 022 ), the effect of the fiber is minimized. 
Here, all of the envelopes are nearly symmetric, and the locations of the transitions between zone A and 
zone B initiation are qualitatively similar between the yield and damage surfaces. 

The effect of clustering is most evident in the initial 012-022 yield and damage envelopes. The predicted 
envelopes appear to be concentric ellipses whose semimajor axes are dictated by the clustering. For 
C, = 0.5, the yield and damage envelopes exhibit similar transitions from zone A to zone B initiation. 
However, for C, = 0. 1 , such transitions are present only in the yield envelope. It is interesting to note that 
while the initial yield and damage points under pure axial shear (which are based on the local fields) are 
reduced by the presence of clustering, the effective axial shear moduli shown in Figure 8 (which are based 
on the two-scale homogenized stiffness tensor) show an increase due to clustering. 

The microstructural details in the case of ca = 0.5 and C = 0.74 and C = 0.24 are shown in Figure 6. 
The corresponding initial yield and damage envelopes are shown in Figure 10. Similarly, for the case of 
ca = 0.74 and C = 0.89 and C = 0.48, the microstructural details are given in Figure 7, and the associated 
initial envelopes are given in Figure 1 1. A comparison of Figure 9, Figure 10, and Figure 1 1, in which the 
volume fraction of zone A, ca, increases from 0.26 to 0.5 to 0.74, reveals that, while all representations of 
clustering provide smaller yield envelopes compared to the uniform case, the ordering of high C, and low C 
envelopes switches. In Figure 9, the higher C, envelopes are smallest, in Figure 10 the higher and lower C 
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Initial Damage Envelopes 
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Figure 9. — Initial yield and damage envelopes for ca = 0.26 associated with Figure 5. The symbols on certain 
envelopes indicate yield/damage initiation within zone A, whereas, when no symbol is present, failure initiated in 
zone B. 


NASA/TM— 2015-218838 


17 


Initial Yield Envelopes 



CT 22/^m 


Initial Damage Envelopes 




jo 

cn 

PO 

b 


— zeta = 0.74 

uniform 3 

— zeta = 0.24 ^ 

// n 



if 

1 -2 I 1 

UdK 2 4 

-2 


-3 


-A- 



<7 22/ <7 cr-n 




1 5 


Figure 1 0. — Initial yield and damage envelopes for ca = 0.5 associated with Figure 6. The symbols on certain 
envelopes indicate yield/damage initiation within zone A, whereas, when no symbol is present, failure initiated in 
zone B. 
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Initial Yield Envelopes 


Initial Damage Envelopes 
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Figure 1 1 . — Initial yield and damage envelopes for ca = 0.74 associated with Figure 7. The symbols on certain 
envelopes indicate yield/damage initiation within zone A, whereas, when no symbol is present, failure initiated in 
zone B. Note that in the 011-022 envelopes for C, = 0.48, under pure on loading, yield and damage initiated in zone B 
(these two points on these envelopes have no symbol). 
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envelopes are very close, and in Figure 11, the lower £ envelopes are smallest. This observation is 
reasonable since the smallest envelopes are always associated with the higher local fiber volume fraction 
being located in the zone with the lower volume fraction within the global RUC. This situation represents 
the greater degree of fiber clustering and leads to higher stress concentrations in the composite, and thus 
contracted initial envelopes. It should be noted that, as before, the transitions of yield and damage 
initiation from zones A to B are denoted in Figure 10 and Figure 11. 

To examine the effect of c A on the initial yield and damage envelopes in a consistent way, the actual C 
value associated with each of the discretizations shown in Figure 4 can be calculated. This is done by 
rearranging Equation (27), 


C act=cj — (30) 

where cj , the average volume fraction in zone A, is calculated for each Ca value by simply averaging the 
fiber volume fractions of the grid squares in zone A. For example, for ca = 0.26, cj is equal to the 
average of the shaded entries in Table 1. For ca = 0.26, C, a ct = 0.35 (cj = 0.62 and c j = 0.40); for 
CA = 0.5, t,ac,= 0.61 (cj = 0.56 and cj = 0.36); for c A = 0.74, £** = 0.83 (cf = 0.51 and cj = 0.30). 

The effect of c A on the predicted composite initial yield and damage envelopes of the carbon/epoxy 
composite, where C, = C, ac t, is shown in Figure 12. Clearly, the variation among the envelopes is less than 
that shown in Figure 9 to Figure 11. This is expected as, for all three cases of c A = 0.26, 0.5, and 0.74, the 
Cad value is between the values of C used in Figure 9 to Figure 11. However, this also indicates that the 
choice of ca used to discretize the microstructure is not critical. 

In order to illustrate the effect of material system on the impact of fiber clustering, the initial yield 
and damage envelopes of an E-glass/epoxy composite have also been examined. The considered isotropic 
elastic properties of the E-glass fibers are: E = 73 GPa and v =0.22, while the same epoxy matrix is 
considered. The mismatch in transverse Y oung's modulus between the fiber and matrix is known to have 
a significant effect on the local fields within composites. For the carbon/epoxy composite considered 
above, this mismatch is 12.4/5.2 = 2.4, whereas for the glass/epoxy composite, it is 73/5.2 = 14. Note that 
the identical carbon/epoxy microstructure (Figure 3) and idealizations into zones (Figure 4) were used. 

Figure 1 3 shows the effect of ca on the predicted composite initial yield and damage envelopes of the 
E-glass/epoxy composite, where C, = C, a ct. Compared to the carbon/epoxy composite (Figure 12), it is clear 
that the effect of clustering on the initial yield and damage envelopes is much greater in the case of the 
E-glass/epoxy composite. The E-glass/epoxy’s greater mismatch in the fiber-matrix transverse Young’s 
modulus leads to higher stress concentrations due to the presence of clustering, which are manifested in 
the predicted envelopes. Similarly, there are now some noticeable variations in the envelopes based on the 
choice of c A , but, as in the case of the carbon/epoxy composite, this difference is significantly less than 
the difference between the uniform case and any of the clustered cases. 
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Figure 12. — Comparison of initial yield and damage envelopes for the carbon/epoxy composite, for different values 

Of Ca With £ = C,act. 
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Figure 13. — Comparison of initial yield and damage envelopes for the E-glass/epoxy composite, for different values 

Of CA With <; = C,act. 
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Nonlinear Effective Stress-Strain Curves 


To examine the effect of fiber clustering on the nonlinear behavior of the composite, the epoxy matrix 
material response was simulated using the nonlinear elastic Ramberg and Osgood (1943) model. The 
uniaxial form of the constitutive relation in this model is given by, 


( _ \ n 


1 Cl 
+ 

1 Cl 
1° 

G 

E E 

v a 0 j 


(31) 


where E is the Young’s modulus and Go and n are two material parameters that can be determined by 
fitting the material nonlinear stress-strain response. A deviatoric multiaxial generalization is given by 
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ij 


1 + v v 
~ a ij ~ ^ a kk°ij + 



\ n—\ 


u eq 
<*0 


(32) 


where v is the Poisson’s ratio and a eq is the equivalent stress (see Eq. (18)). 

In the implementation of Equation (32) within the multiscale HFGMC model described above, it is 
necessary to cast this equation in its incremental form, 


ASy = Sijki A <3 k i 


(33) 


where Sya are the components of the instantaneous compliance tensor, which can be obtained from 
manipulation of the incremental form of Equation (32). Then, inverting S, the local subcell constitutive 
equation (10), including the Ramberg-Osgood nonlinear effects, is obtained. The remainder of the 
incremental multiscale E1FGMC theory remains the same as presented above. 

The nonlinearity of the epoxy material has been modeled herein by choosing two sets of parameters 
with different degrees of nonlinearity: Go = 158 MPa, n = 4 (from Aboudi, 1991) and Go = 80 MPa, n = 15 
(fictional values, resulting in greater nonlinearity). The resulting uniaxial stress-strain curves in simple 
tension are shown in Figure 14. 

The resulting effective transverse stress-strain curves of the composite at the global scale are shown 
in Figure 15 for the three values of ca with varying C and for both sets of nonlinear parameters. It can be 
readily observed that the effect of clustering on the C/epoxy composite is greater when more matrix 
nonlinearity is present. It would be expected that this effect would continue to increase with decreasing 
matrix instantaneous stiffness. In addition, by comparing the three plots in Figure 15, it can be seen that 
for Ca = 0.26, the curve associated with higher C has the greater deviation from the uniform case, whereas 
for Ca = 0.74, the curve associated with lower C has the greater deviation from the uniform case. For 
intermediate value of Ca = 0.5, the curves associates with both values of C are very similar. As in the case 
of the initial envelopes, the greatest deviation from the uniform results is always associated with the 
higher local fiber volume fraction being located in the zone with the lower volume fraction within the 
global RUC. 

Once again, to illustrate the effect of the contrast between the properties of the fibers and the matrix 
on the impact of fiber clustering, the stress-strain response of the E-glass/epoxy composite has been 
examined. The more nonlinear epoxy Ramberg-Osgood parameters. Go = 80 MPa, n = 15, were employed 
in the following simulations. 
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Figure 14. — Uniaxial stress-strain curves of the nonlinear epoxy in 
simple tension, as represented by the Ramberg-Osgood model with 
£ = 5.2 GPa, v= 0.35 and the two sets of nonlinear parameters 
indicated. 



Figure 15. — Effective transverse response of the C/epoxy composite using two sets on nonlinear parameters for 
the matrix, (a) ca = 0.26, (a) ca = 0.5, (a) ca = 0.74. 
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Figure 16. — Effective transverse response of the E-glass/epoxy composite compared with the C/epoxy 
composite, (a) ca = 0.26, (a) ca = 0.5, (a) ca = 0.74. 


The results of the transverse tensile simulations of the E-glass/epoxy composite are compared to the 
C/epoxy composite in Figure 16. Obviously, because the E-glass fiber has a higher transverse stiffness, 
the E-glass/epoxy composite also has a higher effective transverse stiffness. In addition, the higher 
transverse mismatch in fiber-matrix properties results in a much greater influence of fiber clustering on 
the effective transverse response. Otherwise, the trend with respect to the parameter C, is the same in the 
E-glass/epoxy composite as in the C/epoxy system. 

Alternative Microstructure 

To investigate the influence of the assumed global microstructure, an alternative global 
microstructure, taken from a different specimen of the same type of composite, was considered. The 
micrograph is shown in Figure 17. As before, a square grid was placed over the micrograph and the local 
fiber volume fraction within each square was determined via image pixel analysis. Because the aspect 
ratio of the present micrograph shown in Figure 1 7 is different than the previous micrograph shown in 
Figure 3, a grid of 10 x 24 squares was required. The overall fiber volume fraction, c/, determined via 
image pixel analysis of the alternative micrograph is 0.452. A comparison of Figure 5 and Figure 18 
indicates that the two global RUCs are quite different, with the latter having the majority of zone A 
located in the right half of the RUC. 
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Figure 17. — Alternative micrograph from a 0° ply in an IM7/MTM-45 
carbon/epoxy composite. The overall fiber volume fraction, based on 
image pixel analysis, is 0.452. 
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Figure 18. — Multiscale material distribution for ca = 0.25 with C, = 0.5, for the alternative composite 
micrograph shown in Figure 17. 


In order to compare to the previous results, a ca = 0.25 was obtained using a fiber volume fraction cutoff of 
0.55. Note that it was not possible to match the previous ca = 0.26 exactly. The resulting global scale RUC is 
shown in Figure 18, along with the local RUCs and local fiber volume fractions corresponding to C, = 0.5. 

A comparison of the initial yield and damage envelopes of the original microstructure with the present 
alternative microstructure is shown in Figure 19. These comparisons are shown for ca = 0.25 or 0.26 and 
C, = 0.5. As can be seen, although the micrographs used to generate the global RUCs are quite different, the 
predicted initial envelopes under normal loadings are nearly identical. However, for the axial shear case, the 
envelopes show variation under significant shear loading. Thus because of this variation in the shear 
responses, it is not possible to claim that the clustering is fully characterized by the parameters ca and C 
Rather, the chosen global scale microstructure can also have an effect. It should be noted, however, that, 
since only two particular microstructures were considered, the sample size is in no way statistically 
representative. Similarly, no assessment has been made of the statistical relevance of two specific 
microstructures considered in terms of their size and number of fibers included. 
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A comparison between the global transverse stress-strain curves obtained from the original 
microstructure and the present alternative microstructure is shown in Figure 20, in which ca = 0.25 or 
0.26 and C, = 0.5. It can be observed that the two responses are again quite similar. 
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Figure 19. — Comparison of initial yield and damage envelopes for ca = 0.26 and £ = 0.5 from the original (Figure 5) 
and alternative (Figure 18) microstructures. 
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Figure 20. — Comparison of the effective transverse response of the 
C/epoxy composite for the original and alternative microstructures 
with ca = 0.26 and C, = 0.5. The nonlinear epoxy parameters, 
ao = 80 MPa, r? = 15, were employed. 

Conclusions 

A multiscale micromechanical approach has been developed for the prediction of the behavior of 
composites exhibiting fiber clustering. The method is able to predict not only effective composite 
properties, but also initial yield and damage envelopes, as well as the nonlinear stress-strain response of 
the composite, due to the availability of local fields. This is possible because this approach is based on the 
establishment of the concentration tensors on the local and global scales. Fiber clustering is characterized 
by two nondimensional parameters, which measure the extent and density of the fiber clusters, along with 
a global RUC determined from actual composite micrographs. This methodology was applied to 
carbon/epoxy and E-glass/epoxy composites. It has been shown that fiber clustering results in a decrease 
in the sizes of the initial yield and damage surfaces as compared to the uniformly distributed case. The 
effect is greater in the case of E-glass/epoxy compared to the baseline carbon/epoxy due to greater 
mismatch in the fiber-matrix transverse Young’s modulus mismatch in the case of the E-glass/epoxy. In 
addition, transverse stress-strain simulations indicated that the effect of fiber clustering increases with 
increasing fiber-matrix property mismatch and with increasing nonlinearity of the matrix. 

The present method, which can consider three dimensional composite microstructures, could also be 
applied to the simulation of nanocomposites, whose low volume fractions tend to result in significant 
amounts of clustering. Ceramic matrix composites, which typically have brittle matrices and lower fiber 
volume fractions than polymer matrix composites, would also be expected to exhibit significant effects of 
fiber clustering. Further, since FIFGMC has been formulated to include inelasticity, fiber-matrix 
debonding (Aboudi et al., 2013), fiber misalignment (Bednarcyk et ah, 2014), and progressive damage 
(Bednarcyk et ah, 2010; Pineda et ah, 2013), these effects can be incorporated at the local scale in the 
present multiscale implementation to determine the impact of their combination with fiber clustering on 
the composite behavior. 
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